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ABSTRACT 


I* 


NUMESH  is  a computer  program  for  the  numerical 
generation  of  boundary -fit ted  coordinate  systems  for 
two-dimensional  regions.  A numerical  transformation 
maps  a doubly -connected  region  bounded  bv  arbitrary 
curves  onto  a rectangular  computational  field  with  a 
square  mesh.  Numerical  solutions  of  partial  differential 
equations  may  be  obtained  from  the  computational  field 
without  interpolation  regardless  of  boundary  shape.  The 
mathematical  procedure  and  the  use  of  the  program  arc 
described.  Two  types  of  meshes  can  be  produced  which 
are  useful  for  solving  problems  such  as  flow  past  a 
submerged  body  in  a channel  or  under  a free  surface. 

The  use  of  interactive  graphics  permits  a mesh  to  be 
generated,  viewed,  and  regenerated  with  slight  alterations, 
until  a suitable  mesh  definition  is  obtained.  Meshes 
are  shown  for  several  configurations  involving  circular 
and  arbitrarily -shaped  bodies. 


INTRODUCTION 

The  value  of  finite-difference  schemes  for  solving  partial 
differential  equations  has  been  demonstrated  by  their  wide-spread  use  in 
many  areas  of  applied  mathematics,  particularly  in  the  field  of  fluid 
dynamics.  The  usual  problem  is  one  in  which  a partial  differential 
equation  or  system  of  partial  differential  equations  is  to  be  solved  on 
some  bounded  region.  Many  of  the  available  schemes  are  limited  by 
complicated  geometries  which  give  rise  to  difficulties  stemming  from 
inaccurate  numerical  representation  of  boundary  conditions.  The  boundary 
conditions  are  usually  best  represented  when  the  boundaries  themselves 
are  coordinate  lines.  As  a result,  much  use  has  been  made  of  "natural" 
coordinate  systems  such  as  cylindrical  and  spherical  coordinates, although 
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the  practical  usefulness  of  such  systems  is  restricted  to  a few  specialized 
cases.  A general  method  was  needed  which  would  give  a boundary-fitted 
curvilinear  coordinate  system  for  any  specific  geometric  configuration, 
such  as  that  of  an  arbitrary  body  shape  under  a free  surface. 

This  report  describes  NUMESH,  a computer  program  for  the  numerical 
generation  of  such  a coordinate  system  for  an  arbitrary  two-dimensional 
geometry.  NUMESH  uses  a procedure  developed  by  Thompson*  which  leads  to  a 
numerical  transformation  mapping  a set  of  grid  points  in  the  physical 
region  under  consideration  to  a square  mesh  so  that  each  boundary  is 
coincident  with  a grid  line  or  portion  of  a grid  line.  The  numerical 
solution  of  a partial  differential  equation  in  the  physical  region  may 
then  be  calculated  on  the  square  mesh  without  interpolation,  regardless  of 
body  shape  or  mesh  spacing. 

In  Thompson's  early  work,1  he  describes  a mapping  in  which  a body  in 
the  interior  of  a physical  region  is  mapped  to  one  side  of  the  rectangular 
transformed  plane  and  the  outer  boundary  is  mapped  to  the  opposite  side. 

The  remaining  two  sides  of  the  calculating  region  are  the  images,  under 
the  transformation,  of  a cut  running  from  the  outer  boundary  to  the  body. 
The  type  of  coordinate  system  thus  created  has  many  of  the  characteristics 
of  polar  coordinates  and  appears  to  be  suitable  for  problems  in  which 
there  is  only  one  important  boundary,  such  as  fluid  flow  about  a body  in 
an  infinite  fluid. 

2 

Thompson  also  presents  a transformation  in  which  the  body  is  mapped 
to  a slit  within  the  rectangle  and  the  outer  boundary  of  the  physical 
region  is  mapped  to  the  sides  of  the  rectangle.  This  type  of  coordinate 
system  is  more  rectangular  than  the  one  first  described  and  seems  suited 
to  problems  in  which  there  is  more  than  one  important  boundary,  e.g.,  flow 


Thompson,  J.F.,  F.C.  Thames  and  C.W.  Mastin,  "Automatic  Numerical 
Generation  of  Body-Fitted  Curvilinear  Coordinate  System  for  Field  Con- 
taining Any  Number  of  Bodies,"  Journal  of  Comp.  Physics,  vol.  15  (1974), 
pp.  299-519. 

^ Thompson,  J.F.,  F.C.  Thames,  C.W.  Mastin  and  S.P.  Shanks,  "Use  of  Numer- 
ically Generated  Body-Fitted  Coordinate  Systems  for  Solution  of  the  Navier- 
Stokes  Equations,"  Proceedings  of  AIAA  2nd  Comp.  Fluid  Dynamics  Conf . , 
Hartford,  Conn.  (1975). 
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about  a body  in  a channel.  NUMESH  can  produce  this  type  of  mesh,  known 
here  as  Type  1. 

For  the  problem  of  flow  about  a hydrofoil  beneath  a free  surface, 
Thompson  suggests3  another  configuration  in  which  the  transformed  region 
is  made  up  of  two  adjacent  rectangles.  Here  the  body,  the  outer  boundary, 
the  free-surface,  and  the  cut  running  from  the  free  surface  to  the  body 
are  all  mapped  to  the  exterior  of  the  transformed  region.  The  resulting 
mesh  is  used  by  Thompson  for  the  hydrofoil  problem  in  a fluid  of  infinite 
depth. 

NUMESH  can  produce  a new  type  of  mesh,  referred  to  as  Type  2,  which 
has  some  characteristics  of  all  three  of  the  above  configurations.  The 
Type  2 coordinate  system,  described  in  detail  in  this  report,  was  designed 
in  part  for  the  problem  of  an  arbitrary  body  under  a free  surface  in  a 
fluid  of  finite  depth.  Because  it  is  polar  near  the  inner  boundary  and 
rectangular  near  the  outer  one,  the  Type  2 system  allows  the  user  to 
refine  the  mesh  in  the  vicinity  of  the  body  with  little  change  elsewhere 
in  the  field. 


DESCRIPTION  OF  THE  METHOD 


The  brief  discussion  of  the  mathematical  formulation  of  the  method, 
presented  here,  and  the  notation,  follow  those  given  by  Thompson* 3 and 
provide  background  for  the  development  of  the  finite  difference  scheme 
used  in  NUMESH. 

We  wish  to  transform  a two-dimensional , doubly-connected  region  of 
arbitrary  shape  in  the  physical  plane  into  a rectangular  region.  The 
transformation  functions  from  the  physical  plane  (x,y)  to  the  transformed 
plane  (C,n)  are  generated  by  solving  an  elliptic  system  with  boundary 


Thompson,  J.F.,  F.C.  Thames,  J.K.  Hodge,  S.P.  Shanks,  R.N.  Reddy,  and 
C.W.  Mastin,  "Solutions  of  the  Navier-Stokes  Equations  in  Various  Flow 
Regimes  on  Fields  Containing  Any  Number  of  Arbitrary  Bodies  Using  Boundary - 
Fitted  Coordinate  Systems,"  V.  International  Conf.  on  Numerical  Methods  in 
Fluid  Dynamics,  Enschede,  The  Netherlands  (1976). 


conditions  such  that  one  of  the  transformed  coordinates  is  constant  on 
each  of  the  physical  boundaries.  A convenient  choice  for  the  elliptic 
generating  system  is  the  Poisson  equation  because  the  inhomogeneous  terms 
allow  for  some  control  over  the  generated  coordinate  system.  Let  the 
generating  system  be 


Sex  * Sr  = Ci  exp  (‘  ^ (?-^i)2+(n-ni)2  ) = P(C,n) 
>xx  ♦ Hyy  ♦ r Di  exp  |-  / (t-ti)2+(n-ni)2  } h QU,n) 


(1) 


with 

or 


£ = const, 
n = const. 


(2) 


on  each  boundary. 

Since  all  computations  are  to  be  done  on  the  square  mesh  in  the 
transformed  plane,  the  dependent  and  independent  variables  are  interchanged 
giving 


»*«-  28  YXnn  ’ ‘ j2(PxC  * 9V 

«>'«-  2B>v  y>\,„  • • * V 


where 


a 

6 

Y 


x2  ♦ y 


" Vn  + Vn 
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J 


Vn 


- x yr 
n'C 


The  boundary  conditions  are  the  physical  coordinates  of  mesh  points  on  the 
boundaries. 

The  ability  to  alter  the  spacing  of  the  grid  lines  in  the  physical 
plane  is  necessary  to  improve  the  accuracy  of  the  numerical  solution  of  the 
partial  differential  equation  of  ultimate  interest.  Since  the  physical 
coordinates  of  the  points  at  which  the  grid  lines  intersect  the  boundaries 
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are  input  to  the  procedure,  some  control  of  spacing  near  the  boundaries 
can  be  achieved  by  adjusting  this  input.  Further  control  of  the  mesh 
configuration  in  the  field  is  achieved  by  varying  the  F and  Q functions  in 
the  generating  system.  The  f;  and  n lines  may  be  attracted  to  or  repelled 
from  the  specified  points  in  the  field  through  proper  choice  of 

the  and  IT. 

The  approximation  of  Equations  (3)  and  (4)  using  second-order, 
central  differences  for  all  derivatives  yields 


x.  . 
1.3 
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Hquations  (5)  and  (6)  may  then  be  solved  on  the  rectangular  trans- 
formed field.  NUMESH  uses  an  accelerated  Dauss -Seidel  iteration  scheme 
to  obtain  this  solution,  although  other  methods  may  work  as  well.  R,  the 
relaxation  factor  used  to  speed  the  convergence,  must  be  chosen  between 
0 and  2.  The  computation  is  said  to  have  converged  after  the  kth  iteration 


when 


and 


k 

X. 

-x*'1 

< e. 

k 

x.  . 

i.J 

1 

i.J 

k k-1 

y.  . - y.  . 

‘ C1  | 

k 

yiJ 

for  all  i,j  such  that 


and 


> e. 


> e. 


(8) 


There  is  an  optimum  value  of  R,  R^,  in  the  range  0 to  2,  which  causes 
the  scheme  to  converge  in  the  minimum  number  of  iterations.  Experimentation 
with  various  meshes  has  shown  that  this  value  is  often  about  1.6.  A 
choice  of  R greater  than  may  cause  the  iteration  scheme  to  diverge.  The 
iteration  procedure  is  halted  after  the  computation  has  converged  or 
after  a maximum  of  50  iterations.  If  the  convergence  criterion  has  not 
been  met  after  50  iterations,  the  calculation  may  be  converging  slowly  or 
may  be  diverging.  The  printed  output,  described  in  the  following  section, 
should  be  an  aid  in  determining  which  is  the  case.  If  the  computation 
appears  to  be  slowly  converging,  the  relaxation  factor  could  be  increased. 
Additional  iterations  can  be  done  by  using  the  current  results  as  input  for 
another  program  run.  If  the  calculation  appears  to  be  diverging,  all  user- 
supplied  input,  especially  the  relaxation  factor,  should  be  examined  to 
determine  the  cause.  Available  interactive  graphics  routines  are 
particularly  helpful  in  locating  errors  in  the  boundary  input  data. 

As  stated  earlier,  NUMESH  generates  two  types  of  boundary -fit ted 
meshes  when  the  physical  region  is  doubly -connected,  being  bounded  by  two 
arbitrary  non- intersecting  curves.  The  exterior  curve  is  mapped  to  the 
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perimeter  of  the  computation  region  and  the  interior  curve,  hereafter 
referred  to  as  a body,  as  in  fluid  flow  problems,  is  mapped  to  a slit 
within  the  region. 

Type  1 : All  grid  lines  connect  either  two  points  on  the  outer 

boundary  or  a point  on  the  outer  boundary  and  a point  on  the  body.  See 
Figure  1. 

Type  2:  A predetermined  number  of  grid  lines  encircle  the  body  and 

close  on  themselves.  Special  computer  code  is  needed  to  calculate  grid 
points  near  the  body  for  this  type  of  mesh.  See  Figure  2. 

Both  types  of  meshe;  are  suitable  when  accuracy  is  needed  near  both 
the  exterior  and  interior  ncundaries  such  as  in  problems  involving  the 
flow  of  a fluid  past  a submerged  body  in  a channel  or  under  a free 
surface.  A Type  2 mesh  seems  to  be  useful  when  high  accuracy  is 
particularly  important  near  the  body. 

Examples  of  the  Type  1 mesh  are  given  by  Figures  3,  4,  5,  and  6. 

The  Type  2 mesh  is  illustrated  by  Figures  7,  8,  and  9.  The  input  parameters 
discussed  in  the  following  section,  which  were  used  to  generate  these 
examples,  are  summarized  in  Table  1. 
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Figure  3 Type  1 Mesh  About  A Circular  Cylinder  — No  Attraction 


Figure  4 Type  1 Mesh  About  A Circular  Cylinder  - With  Attraction 


— 

Figure  5 Type  1 Mesh  About  An  Arbitrary  Body 
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Figure  6 Type  1 Mesh  Developed  For  Flow  About  A Body  In  A Constricted  Channel 


Figure  7 Type  2 Mesh  About  A Circular  Cylinder  - No  Attraction 


Figure  9 Type  2 Mesh  About  An  Arbitrary  Body 


I 


TABLE  1 - INPUT  AND  CONTROL  PARAMETERS 


Parameter 

Fig.  3 

Fig.  4 

Fig  6 

Fig.  6 

Fig  7 

Fig  8 

Fig  9 

MAXI 

22 

22 

22 

20 

. 

48 

48 

48 

MAXJ 

21 

21 

21 

61 

41 

41 

41 

ISLT 

11 

11 

11 

10 

24 

24 

24 

JLT 

8 

8 

8 

27 

18 

18 

18 

JRT 

14 

14 

14 

41 

24 

24 

24 

El 

01 

01 

01 

01 

01 

01 

01 

E2 

001 

001 

001 

001 

001 

001 

001 

R 

1 6 

16 

16 

16 

16 

1.6 

16 

PLT 

00 

200 

20  0 

25  0 

00 

2 * 103 

2 x 103 

PRT 

00 

20  0 

30  0 

25  0 

00 

2 x 103 

2 x 103 

QLT 

00 

20  0 

200 

25  0 

00 

2 x 103 

2 x 103 

00 

1 

200 

1 

300 

1 

25  0 
1 

0.0 

2 

2 x 103 
2 

2 x 103 
2 

INTL 

0 

1 

0 

0 

0 

1 

0 

NCRIB 

- 

- 

- 

- 

3 

3 

3 

6400  CPU 
TIME 

5 tec. 

4 sec 

7 sec 

25  sec. 

54  sec 

34  sec 

61  sec 
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INPUT 


USE  OF  THE  PROGRAM 


All  input  to  the  program  is  given  in  the  form  of  80-column  card 
images  with  variables  in  free  format. 

The  variables  MAXI,  MAXJ,  ISLT,  JLT,  JRT  are  specified  on  the  first 
data  card.  See  Figures  1 and  2. 

Variable  Description 

MXI  Number  of  mesh  lines  in  n-direction  (SO  max.). 

MAXJ  Number  of  mesh  lines  in  s -direction  (50  max.). 

ISLT  n line  to  which  upper  portion  of  body  is  mapped. 

JLT,  JRT  Extremities  of  body  on  £ axis. 

The  variables  El,  E2,  R are  specified  on  the  next  data  card. 

Variable  Description 

El  Maximum  percent  change  in  values  allowed  for 

convergence,  Ej  in  Equation  (7). 

E2  Minimum  absolute  value  of  x and  included  in 

convergence  test,  Fquati°n  (8). 

R Relaxation  factor. 

The  variables  PLT,  PRT,  QLT,  QRT  are  specified  on  the  succeeding  data 


card. 


Variable 


Description 

PLT  Parameter  controlling  S-line  attraction  to 

points  (JLT, ISLT)  and  (JLT.ISLT+l) . 

PRT  Parameter  controlling  s-line  attraction  to 

points  (JRT, ISLT)  and  (JRT.ISLT+l) . 

QLT  Parameter  controlling  n-line  attraction  to 

points  (JLT, ISLT)  and  (JLT,ISLT-*-l) . 

QRT  Parameter  controlling  n-line  attraction  to 

points  (JRT, ISLT)  and  (JRT.ISLT+l) . 

The  variable  MESH,  control  variable  for  the  type  of  mesh  to  be 
generated,  is  specified  on  the  next  data  card.  Setting  MESH=1  results  in 
a Type  1 mesh;  MESH-2  produces  a Type  2 mesh.  If  MESH=2,  the  variable 
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NCRIB,  the  number  of  grid  lines  completely  encircling  the  body,  is  given 
on  the  card  following.  The  NCRIB  data  card  is  omitted  if  MESH*1 . 

The  control  variable  for  mesh  initialization,  INTL,  is  specified  on 
the  succeeding  data  card.  If  INTL=0,  boundary  input  data  follows 
iimediately.  If  INTL><0,  the  mesh  is  initialized  by  reading  the  previously 
generated  mesh  coordinates  from  input  file  device  TAPE44.  The  boundary 
input  data  cards,  if  needed,  specify  the  x-  and  y-coordinates  of  the 
boundaries.  Each  card  has  an  x value  followed  b>  a y value;  the 
coordinates  are  read  in  the  sequence  given  in  the  following  table. 

(Refer  to  Figure  1 for  Type  1 mesh  and  Figure  2 for  Type  2 mesh.) 


Type  1 (Figure  1) 

1)  A to  B 

2)  D to  C 

3)  A to  D 

4)  B to  C 

5)  E to  F (upper) 

6)  E to  F (lower) 


Type  2 (Figure  5) 

1)  A to  B 

2)  E to  D 

3)  A to  E 

4)  B to  D 

5)  K to  L 

6)  N to  M 


The  following  table  summarizes  the  user-supplied  input  needed  to 
generate  each  type  of  mesh: 

Type  1: 

Card  No.  Variables 

1 MAXI,  MAXJ,  ISLT , JLT,  JRT 

2 El,  E2,  R 

3 PLT,  PRT,  QLT,  QRT 


Type  2: 


boundary  input  data  (if  needed) 

MAXI,  MAXJ,  ISLT,  JLT,  JRT 

El,  E2 , R 

PLT,  PRT,  QLT,  QRT 

MESH 

NCRIB 

INTL 

boundary  input  data  (if  needed) 
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OUTPUT 


The  printed  output  consists  of  all  the  input  parameters,  the  number 
of  iterations  performed,  and  ERROR.  ERROR  is  the  maximum  percent  change 
in  the  x-  and  y-values  that  occurred  during  the  last  iteration.  This 
output  may  be  used  to  check  the  accuracy  of  the  input  data  and  the 
convergence  of  the  iteration  scheme.  The  arrays  X,  Y,  SI,  TA  are  written 
in  free  format  to  the  output  file  device  TAPE33  in  a form  suitable  for  use 
as  initialization  input  to  a subsequent  NUMESH  run.  Initialization  in 
this  way  usually  results  in  fewer  iterations  when  a mesh  is  desired  which 
is  little  changed  from  the  previous  one. 

Variables  Description 

X,  Y Physical  coordinates  of  mesh  points. 

SI,  TA  Parameters  needed  for  potential  flow  solution 

on  generated  mesh. 

I,  J Transformed  coordinates  of  the  mesh  points. 

CONTROL  CARDS 


It  is  usually  desirable  to  save  and  catalog  the  output  from  NUMESH 

for  future  use.  The  control  cards  and  deck  structure  needed  for 

execution  on  the  CDC  6700  are  as  follows: 

J0B  card 
CHARGE  card 
FTN. 

REQUEST ,TAPE33 ,*PF. 

+ ATTACH, TAPE4 4, [previous  output  file  name  1 , 1 IM7XXX . 

LG0. 

CATALOG, TAPE33, [output  file  name] , ID=CXXX ,AC= . 

7/8/9 

source  deck 
7/8/9 

data  cards 
6/7/8/9 


f 

This  card  is  omitted  if  mesh  is  not  to  be  initialized  from  a previously 
generated  mesh. 
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GRAPHICS 


In  order  to  obtain  an  accurate  numerical  solution  to  the  partial 
differential  equation  under  consideration,  a suitable  mesh  must  be 
generated.  It  is  therefore  desirable  to  have  a means  by  which  the  output 
from  NUMESH  may  be  displayed  graphically.  Plotting  routines  with  this 
capability  such  as  IMAGE  (documentation  will  be  available  in  the  near 
future)  are  available  in  the  Computation,  Mathematics,  and  Logistics 
Department.  For  details  on  the  use  of  these  programs,  contact  Code  1843. 

COMPUTER  REQUIREMENTS 

NUMESH  is  designed  to  run  on  the  CDC  6700  computer  system  and  requires 
about  75K  octal  words  of  storage.  It  is  difficult  to  estimate  the 
running  time  because  of  the  large  number  of  factors  involved.  The  computer 
time  required  for  NUMESH  depends  on  the  number  of  mesh  points,  the 
attraction  parameters,  the  convergence  criteria,  and  the  geometry  of  the 
problem.  Execution  times  for  the  figures  are  given  in  Table  1.  NUMESH 
compiles  in  under  15  CP  seconds  on  the  CDC  6400  processor. 


PROGRAM  LISTING 

A program  listing  is  given  on  the  following  pages. 
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ooo  oooo  ooo 


PROGRAM  NUNESM  ( INPUT, OUTPUT, TAPFI I, TAPE44, TAPr? .INPUT) 
common  x (22 ,491 ,y(22,49> ,p(22,49> ,0(22. 49) ,maxi ,hayj ,nit , a<  cp 
I EI.E2.  ISLT.R, SI (22,49)  . TA( 22.49)  ,JlT,  JR T . ME SM . NCR I 9 

•••  TMIS  PROGRAM  NUNER ICALL  V GENERATES  I BODY-PITTING  IP  CM 

CALL  INPUT 

IF (NESM.EQ. 1)  CALL  CONPUTl 
IF (MESH. Ed. 2)  CALL  CONPUT 2 
CALL  OUTPUT 
ENO 

SUBROUTINE  INPUT 

COMMON  X(22 ,49) ,Y(22.49> ,P(22,49) ,0(22, 49) .MAXI ,HAXJ ,NIT ,AERP 
1 El,E2,ISLT,R,SK22.49l  ,TA(22,49)  .JLT.JRT.NESM.NCRIO 

•••  THIS  SUBROUTINE  RE  AOS  INPUT  OATA  AND  COHPUTES  TH* 

ATTRACTION  FUNCTION  VALUES 

RE  AO  (5,* ) MAXI.MAXJ.ISLT, JLT.JRT 
NAXJMlxHAXJ-1 
HA  X IN1*MAX I -l 
REACTS, •!  E1.E2.R 

PRINT  102,HAXI,NAX J.ISLT, JLT,JRT,El,E2.R 

102  FORMAT!/*  MAXI**I5,SX*nAXJ»*I5,5X*ISLT**I5,5x* Jl Tx* 

1 IS.5X* JRTx*I5,5X*E1=*F7. 3.5x*E2**F7. 3,5X*Rx  *C5. i/> 

REA0(5.*>  PLT.PRT, OLT.QRT 
PRINT  lOO.PLT.PRT.QLT.QRT 

100  FORHAT(/1X*PLTx*F10.2,5X*oRT**F10.2.SX*OLT=*pi:. 2, SX  *OPT **F  1 £ 
RE*0(5,*(  MESH 

PRINT  105, MESH 
105  FORMAT ( /II X*NESM=* 131 
IF (MESH. EQ. II  CO  TO  10 
RE  AO ( 5 , • 1 MCRI0 
PRINT  200, NCRI 0 
200  FORMAT (/11X*NCRI8=*I3) 

10  RE  AO  (5,  • 1 I NTL 

PRINT  101, INTL  , 

101  FORMAT(/llX*INTL**IS> 

ISLTP1*ISLT*1 
IF(INTL.EQ.0(  CO  TO  20 
OO  9 1*1, MAXI 

OO  9 J*1 ,N A X J 

9 RE  A 0 ( 44  , • ) X(I,JI,T(I.Jt,SI(I,J)  ,TA(I,J1 ,10UM,  JPUH 

PRINT  104 

104  FORMATI/lOX*  MESH  INITALI2E0  FROM  PREVIOUS  RUN*) 

GO  TO  30 

•••  INITIALI2E  GRID 
20  PRINT  103 

103  F ORMAT ( / 10 X • MESH  INITIALI2EO  FROM  INPUT  OATA*) 

IF (NESN, EO, 2)  GO  TO  70 

I SL  TP1* ISL  T ♦ 1 
I SL  T P2*I  SL  T *2 
00  1 J*1 ,MAXJ 

1 REAO  (5,*)  X(t,J),Y(l,J) 

OO  2 J*1 , M A X J 

2 REAO  (5 , •)  X(MAXl,J) ,Y(MAXI,J) 

00  3 I* 1 , I SLT 

3 REAO  (5,*l  X (1,1  ) , Yd  ,1  ) 

X(ISLTB1,1)*X(ISLT,1» 

Y(ISLTPl,l)*Y(ISLT,l) 

OO  4 I*ISL TP2.NAXI 


6 PEPO  T5,»l  XII,  11.  YU.  It 

no  5 i*i.islt 

s repo  i5,M  xh.mpxji,  v<  i, npxji 

XIISLTP1,HPXJI*X(ISLT,HPX  J1 
TTISLTPl.NPXJl *Y  t ISLT.HPXJt 
00  51  l*ISLTP?,HPXI 
51  PE  *0  «5.  •»  XII.HPxji  ,V| I.  npxji 

no  r j*?.hpxjhi 

no  7 i«?,hpxihi 

* II. J»«* lit  J» 

7 vit,J)*vn.ll 
no  5?  j*jlt,jpt 

5?  PEPO  15, ••  XHSLT.JI.7U5LT.JI 

00  60  J*  JL  T , JP  V 

60  PEPO  T5,»»  XIISLTP1.JI.VIISLTP1.JI 

GO  TO  30 

70  ITP.ISLT-XCPIB 
IBH*ISLT»NCPIB*1 

itphi«itp-i 
IBHP1*IBN»1 
no  7i  j*i,hpxj 

71  PEPO  f 5. •»  Xlt.Jl.TIl. J» 

00  7?  JM.M6XJ 

77  PEPO  15. •»  XTHPXt, J1.7|N»xl, Jl 

00  73  IM.JTPHl 

73  PEPO  15. •»  X(t. ll.7H.lt 

00  731  rUTP.IS* 

731  XII,  ll»7ll.  11*0. 

00  737  I-IBWPl.MXXI 
73?  PEPO  15, •»  Kt.lt.VI!  ,1) 

00  76  1*1,1  TPNl 

74  PEPO  «5.»l  X (t.HAXJI  ,7  (T.HAXJI 

00  761  I*|TP,IBM 

761  XU.HPXJ1«T  U,  NPXJt  *0  . 

00  76?  I*IBPPt,NPXI 
74?  PEPO  (5 . ••  XU.HPXJI  ,7(1, NPXJI 

00  S J*?.MPXJH1 

00  8 I*?,NPXINI 

XU,JI*X(1.J> 

8 YU,Jt*V(I,ll 
00  75  J*JLT,JRT 

75  PEPO  (5,  *1  XIISLT,  JI.UISLT.J) 

00  7 6 J=JLT,J«T 

76  PEPO  15,**  XTISLTPl.Jl.TIISLTPl, J1 

C 

C •••  INITIPLI7E  P PHD  0 PPPP7 

C 

30  JL  T 1*  JL  T 
JPT 1 * JRT 
I PTT1*ISLT 
1PTT?«1SLTP1 
ISTPRT*ISLT ♦? 

IE NO* ISL  T- 1 
JSTPPT*? 

JE  NO*  HPX  JH 1 

IFTHESH.EQ.lt  GO  TO  18 
JLT1* JLT-? 

JPT 1 * JPT ♦? 

IPTTlalSLT-NCPIBU 
IPTT?*ISLT  * NCR I B 
tSTPPT*tSLT*NCPIB*? 

I END* ISL  T-NCPI B- 1 
JSTPPT*JLT 
JENO* JPT 

18  00  6 1*1 ,NPXI 

00  6 J* 1 ,NP  X J 
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PCI.JHQCI, J)*0 
CONTINUE 


6 


•••  REAO  •Tmc T ION  PARAMETERS  ANO  COMPUTE  VALUES 

JCTR* IJLT* JRTI /2 
ISITHIMSIT-I 
ISL  TP?* I SL  T *2 

00  it  i*2« i end 
ri*i 

00  It  J* JSTART.JCTP 
PLTlsPLT 

IFCJ.GT.JLT1)  °lti  = -plti 
RJ*  J 

PCI. J)*PCI, J)*PLT1*EXPC -SORT  I CRI-IATTII  ••2MRJ-JLTH  **2>  I 
11  QCI.J>*QCI,J)-QLT»EXPC-S0RTCCRI-IATTl»**2»CRJ-JLTn»»2)> 

00  12  1*2, I E NO 
RI«I 

00  1?  J=JCTR,JENO 
PRT1*RRT 

IFC  J.LT.  JRU)  PR  T 1 *-PR  T 1 
R J*  J 

PC  It  J»  = P <1  . J)  -PRT1»E*°I  -SORT  C C®I  - 1 ATT  1)  ••  2*  I B J - J RT  H I 

1?  QC I,  J) *QCl , J)-ORT*EXP  C-SQRTC  CRI-IATT1)»»2*CR J-JBTn»»2l I 
00  14  1= ISTAPT.NAXIMl 
R I * I 

00  14  J* JSTART.JCTP 

plti*plt 

IFCJ.GT.JLTH  PLT1*-PLT1 
R J*  J 

PCI.J)*PCI,  J>»PLTl»E*®C-SQRTICRI-IATT?»  CBJ-JLT1)  I I 

14  QCI«J)=QCI . J) »QLT» EXP  C -SORT ( C RI  - I A TT2I  • *2*  C B J- JL  Tl  > • *2)  > 

00  15  1= 1ST  APT .MAT  INI 

R I * I 

00  15  J=JCTR,JENO 
PRT1 =PRT 

IFCJ.LT.JRTH  PRT1«-PRT1 
RJ*  J 

P Cl,  J»*PCI , J)-PRT1»E*PC  -SORT  C CR I - 1 ATT  21  ••  2 ♦ C P J-  J RT  11  • *2  ) C 

15  QC  I,  J)*QCI  , J) *QRT»EXP C-SQRTC  CRI-IATT2) **2»  CP  J-  J»  Tl  I • *2>) 

RETURN 

ENO 

SUBROUTINE  C0NPUT1 

01 NE  NS I ON  X0C22. 49) ,VOI 22,49) 

CONNON  X C22 ,49) ,T  122.49)  .PC22.49) , 9C22, 49)  ,MA«I , NHJ ,NIT ,AEPB, 

1 El,  E2.ISLT.P,  SIC  22.491  ,TAC  22,49)  ,JLT.JPT.Nr^M,Nr,PI9 

MACMlxHAXI-1 
MAX  JH1  = HAX  J - 1 
ISLTH1*ISLT-1 
ISL T PI* I SL  T ♦ 1 
ISLTP2*ISL  T »2 
JLTHlrJLT-l 
JRTP1* JRT»1 

•••  THIS  SUBROUTINE  COMPUTES  THE  X ANO  V COOPINATFS  OF  TMF  MfSM  T V°F  j 
NI T*  0 

00  10  1*1, MAXI 
00  10  J* 1, N AX J 
SICl,J)*TACI,J)*0. 

0 CONTINUE 

00  1 M*1 ,50 
NIT*NITM 

•••  SAVE  OLO  X ANO  V VALUES  FOR  CONVERGENCE  CHECK 


24 


00  2 I-l.NAKI 
00  2 J-V.NAXJ 
I0II.JIXII)  J» 
roii.ji-m.ji 

2 CONTINUE 

c 

C •••  CALCULATE  * ANO  » FOB  UPPEB  HALF  OF  PEG  ION 

C 

00  J t>2,ISLTNt 

00  3 J*2 , HA  X JH1 

xx.ixii,  j»ii-xii,j-ii  i/2. 

YX»IX|I,  J*l»-TII,J-m/2. 

YE»  IYII-1.  JI-YIIM.JI  1/2. 

AL»XE*XE*YE*YE 

«1>I*«E*TKME 

CA«U'U»TI«TI 

AJ2«IXX*YE-IE»YXI 

Xll,J)*AL/2./(AL»GAt*  IX  II  , II.J-ll  I -NC/ I AL»GAI • 

1 IXII-1, J*l)  -I  I l-l . J-tl  -XIIM.  J*ll»XII*t,  J-lt  I /<*. 

2 »GA/2  • / I AL*G  A I * I X 1 1 - 1 . Jl  *X  1 1 * 1 « Jl  I 

3 ♦AJ2*PtI, Jl •**/2./(»L*GAI 

A ♦AJ2*QII»JI **E/2./IAL*GAI 

YU,  JI*AL/2./IAL*GAI*  |Y|t,  J*1I*Y  It  ,J-1I  I-0E/|AL*GAI» 

2 »GA/2. /IAL»GAI*I TI I- l, J> »Y  1 1*1, J » * 

3 ♦AJ2»PII,JI »YX/?. /(AL»CA> 

A »AJ2*QtI.JI »YE/2./ IAL*GA» 

xii, jmxoii ,ji  xii , ji-xoii,ji  i 

YII.JMYOII  ,JI  ♦N*tY|I,  JI-YOII.JI  I 
SI tt*JI>AJ?*QII.JI 
TAII»JI*AJ2*PII»JI 

3 CONTINUE 
C 

C •••  CALCULATE  X ANO  T ON  SLIT  AME  AO  OF  OOOY 

c 

1 * I SL  T 

00  A J*  2 « JL  TH1 

XX* I XII. J»ll-Xtt,J-ll  1/2. 

XE* IXII-l, JI-XII »?, Jl I /2. 

yxmyii,  jui-y  ii,  j-ii  »/2. 

YE*IYIt-t.JI-Yll»?,JI I/?. 

AL*XE*XE»YE*YE 

0E«xx»xe»yx*ye 

GA*XX*XX*YX»YX 

AJ2*IXX*YE-XE»YXI»»2 

XI I.J)«AL/2,/l AL»GA»Y IX II, J»ll»x  II, J-l> l-AE/IALYGAI  • 

1 «XII-l,J»ll-XII-l,J-ll-XII»2,JUt»*II*2.J-lll/A. 

2 »GA/2. / IAL»GAI • (X  I I-1,JI+XII»2,JII 

3 ♦AJ2»P|I, Jl »XX/2./ (AL»GA> 

A »AJ2*QII,JI*XE/2./IAL»GAI 

Yl I, J)«AL/2,/IAL»GAI* IY II , JA1IAY  II  ,J-1I I -IE / |AL«  GAI* 

1 IVII-l,J»ll-YII-l,J-il-YII*2,J*ll»YII*2,J-lM/A. 

2 »GA/2./IAL*GAI*(Yt 1-1. JI*Yt I»2.JII 

3 ♦AJ2»P|I,JI»YX/2./IAL»GAI 

A ♦ AJ2*<JII,JI,YE/2./IAL*GAI 

Xlt,JI«XOII,JI*B»IXII,JI-XOII,JI> 

Ylt,JI*YOII,JI»B»IY|I,J|-YOII,Jll 

Sill, Jl ■ AJ2 *01 1, Jl 

TAII, JI«AJ2*PII, Jl 

XI ISL TP1 , Jl * XI  I,  Jl 

Y IISLTP1.JI *Yl I, Jl 

SIIISLTP1, JI*AJ2»QIISLTP1,  J» 

TAHSLTPl,  J)«AJ2»P|ISLTP1,  Jl 
A CONTINUE 

c 


non  » o o o oo 


•••  C«tCUt*TE  * »N0  X ON  St  IT  9EMIN0  900Y 

00  S J* JNTP1 ,H»X JH1 
xx*ixti, j*t>-xii,j-i>>/?. 

XE*IXIt-t,JI-XIl*?,J>>/?. 

YK*IYII,J*ll-YII,J-lll/E. 

ve«iyii-i, ji-yii*?. ji i/?. 

*L*XE»XE»YE»YE 

«i*»*«*T*»rE 

G»*XX»XX*XX»YX 

» j?«  ixx»xe -xe»yxi ••? 

XI  I,  JM»L/?  ./I*l»G»>«  IX  II,  J»ll  »X  II,  J-ll  l-TE  / l*L»G»>» 

i ixi  i-i«  j»n-v  ii-i,  j-n  -x«i  »z.  j*ii»»iw,j-i »» /«• 

Z ♦ &*/? . / <»l*G»> • I XI  I- l, J>  »X I I*?,Jl> 

t ♦ »J?»QII,JI,XC/E.FI*I.«G»> 

XII.  Jt»ftL/Z./m»Gft>«  IX  1 1,  J*ll* XII. J-ll  >-1c/l»l*G»>» 

i «xii-i,j*n-xii-t,j-n-xii*z,j*ii*x(t*?.j-iu/«> 

Z ♦ G»/?./UL»G*l»IVII-l,  JI*XI 

J ♦ 8J2*PII,JI»YX/?./Ut.*G*> 

G ♦»J?»QII,  Jl  *YE/2./ I*l»G*» 

xii.j>*xoii,j>»p»txti,j>-xoii,jn 
XI  I,  JI*Y0I I ,J>  ♦R*IYII,  J1-X0II,  Jl  I 
SI  1 1 . J) * *J?  *Q  1 1 , Jt 
T* (I, J>»* J?*PI I, J> 

XIISLTP1 ,JI xxil, Jl 
YllSlTPl.J)  PXII, Jt 
SI  1 1 St  TP  1 « J)*AJ2*Q( ISLTPl.JI 
T8IISITP1, JI*8J2»P|ISLTP1, Jl 
CONTINUE 

•••  C»lCUl»TE  X »N0  X 9EL  ON  TOOT 

00  G I=ISlTP2»HiXIMl 
00  G J*?.N»XJH1 
XKx|Xlt,J*ll-XtI,J-l>>/?. 

xe«ixii-i, j i -xi i*i ,ji i/e. 
xx»ixn,  J»n-xii,  j-i>  »/2. 

Xf*  IXII-l, JI-XII»l. ji * /2. 

*l*XE»XE»XE»YE 
9E*XX»XE  *x  x»TE 
G»*XX»XX»YX*XX 
»J2»IXX»XE-XE»XX> ••? 

XII,JI*»LF2./UL»G»>*IXIt,J*ll»Xll,J-l> >-9E/I«L»G«>* 
1 IX<J-|,J»ll-X<I-l.J-l!-X!I»t,J*l!»XI!*l,J- 111/4 

Z •G8/?./tftl»G* 1*1X11-1, J)»Xtt*l,Jtt 

J ♦»J?*P(I, J> »XX/?./ I»L»G*> 

fc  ♦»JZ*Qa,JI*XC/?./t«U»G»> 

YII.JI*»L/2./l»L*&»>*IYII, J*ll »XII,J-1> >-1S/I»l*G»>» 
1 IXII-1,J*1> -XII-1. J-U-YIIM, J»ll* YII*1 ,J-1>> /4 

Z ♦G»/2./l»t»G»>»IYI I- 1, J>  »X II»1,JI) 

3 ♦»J?»PII, J> •YX/2./UL»G»> 

G *»J2*QII, Jl •XE/?./I*L»G«> 

XII, JI*XOII,JI ♦*•1X11, Jl-XOfl.JII 
XII, JI«YO!I . J> ♦*•1X11 , JI-XOII.JI I 
Sill. JI*»J2»QII, J> 

T»!I, JI>«J?*PII, Jl 
CONTINUE 

•••  CHECKING  CONVERGENCE  

«€**•!. 

00  T 1*1, NG  XI 
00  T J.JLT.JPT 

IF  1 80S  IX II.JII .IE  • E?l  GO  TO  8 
ENNPMSIlTOIt.JI-TIt,  JII/TII.JII 


IFIERR.GT.AERRI  AERR-ERR 

8 GO  TO  7 
ER»*A8SHX0II.  J>  -X  If.  Jll/X  II,  jll 
IFIERR.GT.AERRI  AERR*ERR 

T CONTINUE 

tFIAERR.LE.Ell  GO  TO  9 

1 CONTINUE 

9 RETURN 

ENO 

SUBROUTINE  COMPUTE 
DIMENSION  XOI22.49) .V0I22.49I 

CONMON  XI22.49l.ri22.%9l,Pt22.49>.9t22.49t , M A X I , MA  X J , Nl  T , AE PP  , 

1 E1.E2.  ISL  T,R,$II22,<»9>  . TAI22.49I  . JLT,  JRT.MESH.NCRIB 

I T®« ISLT -NCRIB 
IBNaISLT*NCRIB»l 
MAXINl*NAXI -1 
MAX JM1.MAX J-l 
ISLTM1*ISLT-1 
ISLTP1.ISLT»1 
ITPNl«ITP-t 
ITPP1MTPM 
I8NP1 * I BN* 1 
I9NH1«IBN-1 
JRTP1*JRT*1 
JLTHlaJlT-1 

C 

C •••  THIS  SUBROUTINE  COMPUTES  THE  X ANO  V COORINX TES  OP  THE  HESH  TVPf  2 
C 

NIT»0 

00  200  I* 1 * MAX  I 
00  200  J*1 « MAX J 
SI  II. JI*TA  (I.Jt-O. 

200  CONTINUE 

00  1 N* 1.50 
NI T»NIT ♦ 1 
C 

C •••  SAXE  OLO  X ANO  T VALUES  FOR  CONVERGENCE  CHECK 
C 

00  2 1*1. MAXI 
00  2 J-l.MAXJ 
xo(I.Jt*x<i,j> 

VOtt.JlaVII.JI 

2 CONTINUE 
C 

c •••  UPOATE  OUNNV  VALUES 
C 

00  10  J»  1. JLTN1 
XIITP, J» *XI IBNP1, J) 

VIITP,JI*VIIBNP1,JI 

10  CONTINUE 

00  11  J*  JR  T PI  * MAX  J 
XIITP. JI*XIIBMP1,JI 
VIITP.JlaVIIBMPl.JI 

11  CONTINUE 
C 

C •••  CALCULATE  X ANO  V FOR  UPPER  HALF  OF  REGION 
C 

00  J Is 2 . 1 T PHI 

00  3 «l* 2.MAXJM1 

XXalXII.  J»1I-XII,J-1ME2. 

XE»IXII-1, JI-XII.l ,J» 1/2. 

VX-IVII,  JR1I-TII.J-1M/2. 

VE-IVIf-l,  JI-VII.l.JM/2. 

AL»XE»XE*VE*VE 
BE*XX*XE ♦VX»VE 
CA*XX*XX.VXaVX 


t m i-i.j»u-« n-t, j-ii -iii»i,j«ii»*ii»i, j-iii/a. 

3 ♦AJ2*PII,JI***/2./IAL*«AI 

% Jl •tl/t./ IAL»GAI 

I-9E  /IAL»GAI* 

1 ivii-i,j*n -v  il-t,  J-n-vii*i»  j*ii*viim,j-i»ma. 

2 «GA/2.  /IAL«GAI»I  VI 1-1,  Jl  »V  I IM.JH 

3 ♦»J?»»II, J» •V*/?./|»L»G»» 

* ♦AJ2»QII,JI*VF/2./IAL»GAI 

*11. JI**OII,J>  »»aIXII ,J»-*0II, Jl  I 

vn,  ji.vot  i,ji*»»ivii,  ji-voii.jii 

Sill, JI*AJ2*0II,JI 

iiii.jiujihii.ji 

3 COntmuE 

c 

c •••  U*OATE  OUHHV  VALUES 
C 

N*0 

00  12  t*ITP,IBH 
H*H*1 

*11,  JLtnil  «*II9HPl-n.  JLTt 
*11,  JATP1I  **IIBHPl-n,  J*tl 
v 1 1, JLTH1I  *V II9HP1 -N, JL T I 
v 1 1,  J»TP1I  *V II9nPl-n,  JPTI 

12  continue 

c 

C •••  CALCULATE  * Ano  V »»OUAD  THE  900V 
C 

00  16  IMTP.IBn 

TFii.ea.iSLT.on.i.eo.iSLTPn  go  to  16 
00  13  J*  JL  T , J9T 

«i*i«ii,  j*n-*  it,  j-m /2. 

*e*  i*ii-i,ji-*ii»i.ji t/2. 

V**IVII,  JHI-VII.J-11  1/2. 

ve*ivii-i,j»-vii*i,jii/2. 

al*xe**e*ve*ve 

BE»XX*KE»VK*TE 

GA=*****»V*»VX 

AJ2*IX**VE-XE*VXI**2 

*11, J»*AL/2./IAL»GAI* l*<t,J«ll»* II,J-1>  l-95/IA(.*G»>» 

1 1*11-1.  JV1I -KII-1. J- II  -*IIM,  JMMXIIM.J-1II/A. 

2 ♦GA/2./IAL*GAI*<* U-1,JI»*II»1,JM 

3 ♦AJ2*PII, Jl ***/2./ IAL »GA) 

A ♦ AJ?,(|II,JI*XF/2./IAL*GAI 

nr. ji*al/?./iai»gai*  it ii, j»ii ♦*  ii ,j-n  i-ie/ ial»g*i» 

1 ivii-i,j»ii-vii-i, J-H -VII *i, j«n»vii»i, J-1II  / A. 

2 »GA/2. / |AL»GAI»IVII-l,JI*V(I»t,JII 

3 »AJ2*P|I,JI»V*/2./|»L»GAI 

A ♦AJ?,0II,JI*VE/2./IAL*GXI 

*11, J1*K0II, Jl »*»I*II, JI-*OII,J» I 
VII, J»*VOII ,JI VII, JI-VOII.JI I 
SIIt,JI*AJ?*OII,JI 
TAII,JI*AJ2*Pll,JI 

13  continue 

16  continue 

c 

c •••  UPDATE  ounnv  VALUES 
C 

00  1A  J* l, JL TM1 
*iien,ji*K  HTPnt,  Jl 
vi I8n. ji  *v I ITPni, j> 
ia  ContmuE 

00  15  J»JPTPl,nA*J 
*il9n,j»*»iiTPni,ji 
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HIBH,J1*Y(ITPN1  ,J1 
19  CONTINUE 
C 

C •••  CALCULATE  I ANO  T BEL  ON  BOOT 

C 

DO  B I*IBNP1,MAXIK1 

00  B J*;,NAIJN1 

n*  mi.  jmi-hi.j-iii/;. 

XF*(I(I-1. Jl-Xlt»i ,J) I/?. 

TI*mi.  4*1 l-Ht,  J-ll  1 /*. 

TE* (TII-1, Jl -T(t*l  .Jl  Id. 

AL*IE*IC*TC»YC 
BE ■ H*IE*T  I*YE 
CA*XX*II»YI*YI 
AJI* (II»YE-IE»YII**; 

X(I.J>*AL/;./IAL»GA>*(Id,JHI»I  II.J-ll  Mt/IAl*6ll* 

i (id-i,  j*ii-hi-i,  j-ii-iiih.j»u»hi«i,j-iii/9 

t •GA/;./IAL«GAI*|I(I-l.  Jl  ♦■  ( Kl.JII 

5 ♦ AJC'Pd.  Jl *X */?./ (AL’CAI 

9 »AJ?*Q(I.JIaIE/;./(AL»CAI 

YII»JI*AL/;./(AL*GAl*(Y(I»J*ll*YII»J-lll-9E/(AL»GAI* 
1 (V(I-l,J«lt-Yd-l.  J- It -HIM.  J«lMV(I»l,J-ltl/<. 

; »GA/;./IAL*GA>*(VII-l,JI*YII»l,JII 

3 ♦Aj;*P(I,Jl*YX/;./(AL*GAI 

* »aj;*q(i,ji*ye/;./ial»gai 
HI.  J 1*1 0(1.  JI«**(HI.JI-X0(I,J>  I 
TII, JI*VO(I ,J| *9*1 Ttl, JI-TOII, Jl I 
SKI.  JI*Aj;*Q(I,  Jl 
TAII, JI*Aj;*B(I, Jl 
B CONTINUE 

C 

C •••  CHECKING  CONVENGENCE 

C 

A£99*8. 

oo  r i*;,haiihi,; 

oo  r J*JLT,J*T 

IF  I ABS IT (I.JII .LE • EEI  GO  TO  A 
ENN*  ABS  (IYOII.JI-Y (I.JII/TII.JII 
IFIEKA.GT.AEAKI  ACNN*  CNN 
a triA8Sfifi.jii.LE.£;i  go  to  r 

ENN*  ABS  MXO(t.  Jl  -1(1,  Jl  l/XII.JI  1 
IF(E**.GT.  AENNI  IE«»*fW 
T CONTINUE 

IFIAEM.LE.E1I  GO  TO  9 
1 CONTINUE 

C 

C •••  UPDATE  DUNNT  VLAUCS 
C 

9 N*0 

00  i;  I*ITP,IBN 
N«N«1 

1(1, JLTH1I*IIIBHPI-N, JLTI 
HI,  JNTP1I  *1(1  BNP1  -N,  JNTI 
T(I,JLTN1I*TIIBHP1-N, jlti 
T(I,JNTP1I*T( IBNP1 -N, JNTI 
IT  CONTINUE 

jltn;*jlt-; 

jhtp;*j*t*; 

00  19  I*ITPP1,IBNN1 

oo  19  j*i,jltn; 

HI.  Jl*0. 
t9  t(i,ji*b. 

00  18  I* ITPP1 , IBNNl 
00  IS  J* JBTPI.NAIJ 
HI.JI-B. 

18  T( I , Jl*0 . 
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RETURN 
END 

SUBROUTINE  OUTPUT 

CONNON  *«!*,%R».TIEE,RRI,RI*2,*9».QI?E.kH.H»*I,>«tJ,mT,»E0®. 

1 El.E2,ISLT,R,sn2?,M»  ,T*IE?.«»RI  ,JLT, JRT.NESN.NCRIB 

•••  TNIS  SURROUTINE  PRINTS  OUTPUT  D»T»  ANO  NRtTES  X ANO  T TO  PILE 

PRINT  t.NIT, AERN.EttE? 

• EORNAT I //*  ITERATIONS«MS,5*»ERROR»»F10.6.SX»E1«*F&. J.5X»E2«»Ff,.  I» 

DO  It  I-l.NAKI 
00  It  JM.NAXJ 

NRITE  (33**1  X«I,J»  ,T«I,  JI.SHI.  J>.T*1I.  J>.I.J 
SI  CONTINUE 

RETURN 
ENO 
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I 


FILMED 


DTNSRDC  ISSUES  THREE  TYPES  OF  REPORTS 

(1)  DTNSRDC  REPORTS,  A FORMAL  SERIES  PUBLISHING  INFORMATION  OF 
PERMANENT  TECHNICAL  VALUE.  DESIGNATED  BY  A SERIAL  REPORT  NUMBER 

(2)  DEPARTMENTAL  REPORTS,  A SEMIFORMAL  SERIES,  RECORDING  INFORMA 
TION  OF  A PRELIMINARY  OR  TEMPORARY  NATURE.  OR  OF  LIMITED  INTEREST  OR 
SIGNIFICANCE.  CARRYING  A DEPARTMENTAL  ALPHANUMERIC  IDENTIFICATION 

(3)  TECHNICAL  MEMORANDA.  AN  INFORMAL  SERIES  USUALLY  INTERNAL 
WORKING  PAPERS  OR  DIRECT  REPORTS  TO  SPONSORS.  NUMBERED  AS  TM  SERIES 
REPORTS.  NOT  FOR  GENERAL  DISTRIBUTION 


